1 Load packages

library(tidyverse)
library(emmeans)
library(glmmTMB)
library(kableExtra)
library(lme4)
library(lmerTest)

2 Load dataset

data <- read_csv("SB12.csv")

2.1 Glimpse dataset

glimpse(data)
## Rows: 160
## Columns: 19
## $ location                 <chr> "Lancaster", "Lancaster", "Lancaster", "Lanc…
## $ year                     <dbl> 2019, 2019, 2019, 2019, 2019, 2019, 2019, 20…
## $ herb                     <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
## $ other                    <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
## $ trade                    <chr> "PRE only", "PRE only", "PRE only", "PRE onl…
## $ plot                     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ trt                      <dbl> 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4,…
## $ rep                      <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3,…
## $ waterhempcontrol_14      <dbl> 85, 85, 95, 35, 99, 99, 99, 99, 91, 95, 85, …
## $ grasscontrol_14          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ waterhempcontrol_28      <dbl> 88, 25, 95, 10, 95, 99, 99, 99, 98, 88, 75, …
## $ grasscontrol_28          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ waterhempcontrol_harvest <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ counts_m2                <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ biomass_gm2              <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ yield_bu                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ phyto_14                 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ dicambadrift_14          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
## $ canopyclosure_14         <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …

2.2 Skim dataset

skimr::skim(data)
Data summary
Name data
Number of rows 160
Number of columns 19
_______________________
Column type frequency:
character 4
numeric 15
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
location 0 1 8 9 0 2 0
herb 0 1 8 46 0 10 0
other 0 1 8 46 0 10 0
trade 0 1 5 14 0 10 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
year 0 1.00 2019.50 0.50 2019.00 2019.00 2019.50 2020.00 2020.00 ▇▁▁▁▇
plot 80 0.50 255.50 112.55 101.00 178.25 255.50 332.75 410.00 ▇▇▁▇▇
trt 0 1.00 5.50 2.88 1.00 3.00 5.50 8.00 10.00 ▇▇▇▇▇
rep 0 1.00 2.50 1.12 1.00 1.75 2.50 3.25 4.00 ▇▇▁▇▇
waterhempcontrol_14 0 1.00 92.71 12.59 35.00 92.00 99.00 99.00 100.00 ▁▁▁▁▇
grasscontrol_14 124 0.22 87.36 10.73 65.00 79.00 89.50 98.00 100.00 ▂▃▃▃▇
waterhempcontrol_28 0 1.00 85.47 19.62 0.00 80.00 94.00 99.00 100.00 ▁▁▁▂▇
grasscontrol_28 120 0.25 90.65 25.88 0.00 100.00 100.00 100.00 100.00 ▁▁▁▁▇
waterhempcontrol_harvest 80 0.50 78.47 28.66 0.00 75.00 88.50 95.00 99.00 ▁▁▁▂▇
counts_m2 80 0.50 3.05 3.55 0.00 1.00 2.00 4.00 18.00 ▇▂▁▁▁
biomass_gm2 80 0.50 17.82 51.05 0.00 0.00 1.92 12.46 382.06 ▇▁▁▁▁
yield_bu 52 0.68 59.21 11.57 30.44 49.23 62.93 68.64 74.80 ▁▃▃▃▇
phyto_14 85 0.47 1.93 2.81 0.00 0.22 1.00 3.00 12.00 ▇▁▁▁▁
dicambadrift_14 120 0.25 2.15 4.45 0.00 0.00 0.00 3.00 25.00 ▇▁▁▁▁
canopyclosure_14 120 0.25 63.38 12.69 20.00 60.00 70.00 70.00 77.00 ▁▁▁▅▇

3 Data visualization

3.1 Waterhemp control at 14 DAT

ggplot(data, aes(x=reorder(other,waterhempcontrol_14), y=waterhempcontrol_14,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

3.2 Waterhemp control at 28 DAT

ggplot(data, aes(x=reorder(other,waterhempcontrol_28), y=waterhempcontrol_28,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

3.3 Waterhemp control at harvest

ggplot(data, aes(x=reorder(other,waterhempcontrol_harvest), y=waterhempcontrol_harvest,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

3.4 Weed biomass

ggplot(data, aes(x=reorder(other,biomass_gm2), y=biomass_gm2,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Weed biomass (g m2)") +
  theme_minimal() +
  theme(legend.position = "none")

3.5 Yield (bu/acre)

ggplot(data, aes(x=reorder(other,yield_bu), y=yield_bu,
       fill=trade, color=trade)) +
  geom_boxplot(color="black") +
  geom_jitter(alpha=0.2) +
  facet_grid(year ~ location) +
  coord_flip() +
  labs(x="", y="Waterhemp control (%)") +
  theme_minimal() +
  theme(legend.position = "none")

# Data wrangling

new_dt <- 
  data %>% 
  rename (herbicide = other) %>% 
  mutate(
    wt_14 = waterhempcontrol_14/100,
    wt_28 = waterhempcontrol_28/100,
    yield_kg = yield_bu * 67.5) %>% 
  mutate(
    wt_14 = 
      case_when(
        waterhempcontrol_14 == 100   ~  0.99,
        TRUE            ~ wt_14),
    wt_28 = case_when(
        waterhempcontrol_28 == 100   ~  0.99,
        waterhempcontrol_28 == 0   ~  0.01,
        TRUE            ~ wt_28)) %>% 
  filter(
    herbicide != "PRE only"
  )

4 Data analysis

4.1 Waterhemp control at 14 DAT

4.1.1 Model

  • Model using herbicide trts, year, and location as fixed effects. Only rep as random effects
model_14 <- glmmTMB(wt_14 ~ herbicide * location * year + (1|rep), beta_family(link = "logit"), data=new_dt)

4.1.2 Anova

glmmTMB:::Anova.glmmTMB(model_14)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: wt_14
##                             Chisq Df Pr(>Chisq)    
## herbicide                 15.4841  8    0.05039 .  
## location                   1.9251  1    0.16530    
## year                    9714.1367  1    < 2e-16 ***
## herbicide:location         6.4935  8    0.59213    
## herbicide:year             1.4430  8    0.99361    
## location:year              0.3624  1    0.54716    
## herbicide:location:year    0.6957  8    0.99954    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Anova shows no interaction of year, herbicide and location. There is only herbicide and year effects

4.1.3 New model

Model using herbicide trt as fixed, and rep, year, and location as random effects

model_14_2 <- glmmTMB(wt_14 ~ herbicide + (1|rep) + (1|location/year), beta_family(link = "logit"), data=new_dt)

4.1.4 New Anova

glmmTMB:::Anova.glmmTMB(model_14_2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: wt_14
##           Chisq Df Pr(>Chisq)
## herbicide 11.21  8     0.1901
  • No effects of herbicide treatments. However, I am looking into herbicide effects.

4.1.5 Least square means

lsmeans_14<- emmeans(model_14_2, ~ herbicide, cont="pairwise", adjust="none", type="response", alpha=0.05)

4.1.6 Compact letter display

plot(lsmeans_14, ~ herbicide, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

CLD(lsmeans_14$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)
herbicide response SE df lower.CL upper.CL .group
3 PRE fb glufosinate + fomesafen + acetochlor 0.9663473 0.0092718 131 0.9423191 0.9805724 a
9 PRE fb glufosinate + fomesafen 0.9644099 0.0096027 131 0.9396867 0.9792228 a
1 PRE fb glufosinate + fomesafen + S-metolachlor 0.9633468 0.0100886 131 0.9372445 0.9788374 a
2 PRE fb glufosinate + pyroxasulfone 0.9513317 0.0122489 131 0.9205233 0.9705792 ab
7 PRE fb glufosinate + imazethapyr 0.9491984 0.0129393 131 0.9165874 0.9694842 ab
8 PRE fb glufosinate + acetochlor 0.9470465 0.0131010 131 0.9142870 0.9677273 ab
6 PRE fb glufosinate 0.9426195 0.0139116 131 0.9080580 0.9646943 ab
4 PRE fb glufosinate + S-metolachlor 0.9408544 0.0142213 131 0.9056144 0.9634678 ab
5 PRE fb glufosinate + dimethenamid-P 0.9296513 0.0162142 131 0.8900147 0.9557144 b

4.2 Waterhemp control at 28 DAT

4.2.1 Model

  • Model using herbicide trts, year, and location as fixed effects. Only rep as random effects
model_28 <- glmmTMB(wt_28 ~ herbicide * location * year + (1|rep), beta_family(link = "logit"), data=new_dt)

4.2.2 Anova

glmmTMB:::Anova.glmmTMB(model_28)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: wt_28
##                            Chisq Df Pr(>Chisq)    
## herbicide                23.1383 10   0.010247 *  
## location                 14.2304  3   0.002608 ** 
## year                    194.8150  1  < 2.2e-16 ***
## herbicide:location       13.3607  8   0.100027    
## herbicide:year            2.3217  8   0.969526    
## location:year             0.0097  1   0.921564    
## herbicide:location:year   0.0139  8   1.000000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Anova shows no interaction of year, herbicide and location. There is only herbicide and year effects

4.2.3 New model

Model using herbicide trt as fixed, and rep, year, and location as random effects

model_28_2 <- glmmTMB(wt_28 ~ herbicide + (1|rep) + (1|location/year), beta_family(link = "logit"), data=new_dt)

4.2.4 New Anova

glmmTMB:::Anova.glmmTMB(model_28_2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: wt_28
##           Chisq Df Pr(>Chisq)    
## herbicide 33.38  8  5.259e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • No effects of herbicide treatments

4.2.5 Least square means

lsmeans_28 <- emmeans(model_28_2, ~ herbicide, cont="pairwise", adjust="none", type="response", alpha=0.05)

4.2.6 Plot

plot(lsmeans_28, ~ herbicide, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

4.2.7 Compact letter display

CLD(lsmeans_28$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)
herbicide response SE df lower.CL upper.CL .group
1 PRE fb glufosinate + fomesafen 0.9335527 0.0219066 131 0.8747880 0.9658156 a
3 PRE fb glufosinate + fomesafen + S-metolachlor 0.9316031 0.0228235 131 0.8702296 0.9651143 a
7 PRE fb glufosinate + fomesafen + acetochlor 0.9297294 0.0230599 131 0.8681042 0.9637635 a
9 PRE fb glufosinate + acetochlor 0.9012060 0.0304182 131 0.8227128 0.9471779 ab
8 PRE fb glufosinate + pyroxasulfone 0.8942952 0.0319956 131 0.8124210 0.9429429 ab
2 PRE fb glufosinate + S-metolachlor 0.8814474 0.0351589 131 0.7925916 0.9353419 ab
5 PRE fb glufosinate + imazethapyr 0.8442439 0.0431294 131 0.7391023 0.9120553 bc
6 PRE fb glufosinate + dimethenamid-P 0.8436929 0.0424061 131 0.7407439 0.9106904 bc
4 PRE fb glufosinate 0.7817276 0.0528757 131 0.6598743 0.8686175 c

4.3 Yield

4.3.1 Yield raw data distribution

library(ggridges)
new_dt %>% 
  ggplot(
    aes(y=herbicide, x=yield_kg, fill=herbicide)) +
      geom_density_ridges(scale=2, show.legend = FALSE)

4.3.2 Data manipulation

yield_dt <- new_dt %>% 
  mutate(
    siteyr = case_when(
      location == "Lancaster" & year == 2020 ~ "Lan20",
      location == "Lancaster" & year == 2019 ~ "Lan19",
      location == "Brooklyn" & year == 2020 ~ "Bro20",
      location == "Brooklyn" & year == 2019 ~ "Bro19",
      TRUE                     ~  "siteyr"
    )
  )

4.3.3 Model

fit <- lmer(yield_kg ~ herbicide * siteyr + (1|rep), data=yield_dt)

4.3.4 Anova

anova(fit)
## Type III Analysis of Variance Table with Satterthwaite's method
##                    Sum Sq  Mean Sq NumDF  DenDF  F value Pr(>F)    
## herbicide          428797    53600     8 67.094   0.5795 0.7912    
## siteyr           35541129 17770564     2 67.950 192.1351 <2e-16 ***
## herbicide:siteyr  1318370    82398    16 67.087   0.8909 0.5820    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.3.5 Least square means

lsmeans_fit <- emmeans(fit, ~ siteyr, cont="pairwise", adjust="none", alpha=0.05)
## NOTE: Results may be misleading due to involvement in interactions

4.3.6 Plot

plot(lsmeans_fit, ~ siteyr, comparisons=TRUE, type="response", alpha=0.05, adjust="none")

4.3.7 Compact letter display

CLD(lsmeans_fit$emmeans, alpha=0.05, Letters=letters, adjust="none", reversed = TRUE) %>%
  kbl() %>%
  kable_classic_2(full_width = F)
siteyr emmean SE df lower.CL upper.CL .group
2 Lan20 4672.755 84.84491 4.867462 4452.856 4892.654 a
1 Bro19 4348.885 94.32107 6.946555 4125.503 4572.267 b
3 Bro20 3311.142 84.84491 4.867462 3091.242 3531.041 c